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We study a linear ramp ol the nearest-neighbor tunneling rate in the Bose-Hubbard model driving 
the system from the Mott insulator state into the superfiuid phase. We employ the truncated Wigner 
approximation to simulate linear quenches of a uniform system in 1,2, and 3 dimensions, and in a 
harmonic trap in 3 dimensions. In all these setups the excitation energy decays like one over third 
root of the quench time. The — | scaling arises from an impulse-adiabatic approximation - a variant 
of the Kibble-Zurek mechanism - describing a crossover from non-adiabatic to adiabatic evolution 
when the system begins to keep pace with the increasing tunneling rate. 

PACS numbers: 03.75.Kk, 03.75.Lm 



INTRODUCTION 



Gaplcss quantum critical points are a serious obsta- 
cle for quantum simulation with ultracold atomic gases 
or ion traps, where one would like to prepare a simple 
ground state of a simple initial Hamiltonian and then 
drive the Hamiltonian adiabatically to an interesting fi- 
nal ground state. This general observation has been re- 
cently substantiated by a more quantitative theory [l], 3 , 
that is by now confirmed by several numerical studies 3[ . 
The theory is a quantum generalization of the classical 
Kibble-Zurek mechanism (KZM) H, The theory pre- 
dicts that density of excitations (or excitation energy) 
decays with (usually a fractional) power of quench rate 
1/tq. Experiments dedicated to the quantum theory 
were made in Refs. [fUl]- In Ref- ultracold spin- 
less bosonic atoms in a three-dimensional (3D) optical 
lattice were driven from the Mott insulator phase to the 
superfiuid phase across a quantum phase transition. The 
transition was non-adiabatic and the excitation energy 
was reported to decay with approximately the third root 
of the transition time. Much in the same spirit, adia- 
baticity of loading atoms into an optical lattice or re- 
leasing them from the lattice confinement [l(| has been 
questioned recently. The 1/3-scaling reported in the ex- 
periment @ coincides with our earlier prediction in Ref. 
fill ] for ID chains. It is the aim of this paper to show 
that the same scaling also holds in 3D. 

In this paper we generalize the theory developed in 
Ref. [ll[ for a ID ring of BEC's to a two- and three- 
dimensional optical lattice in a harmonic trap potential. 
The system is initially prepared in the ground state deep 
in the Mott regime. Since atoms are localized in this 
regime, there is definite number of particles at each site 
that translates into an indefinite phase. The random 
phases at different sites are uncorrelated. After prepa- 
ration of this state, the strength of the lattice potential 
is gradually reduced diminishing potential barriers be- 
tween sites and increasing the nearest-neighbor tunnel- 
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FIG. 1: The transition rate (dJ/dt)/J — J~ Tq is huge at 
early time and small at late time. In contrast, the relaxation 
rate r _1 = J 1 ^ 2 is negligible at early time and large at late 
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time. The two rates are comparable near J 

J the evolution is impulse, i.e., the state does not change and 
remains the initial Mott state with random phases. After J 
the evolution becomes adiabatic. 



ing rate. As the tunneling rate is ramped up towards the 
superfiuid phase, the random initial phases become in- 
creasingly correlated. The range of these correlations as 
well as the excitation energy depend on the rate of this 
transition. We argue that the time evolution with the 
increasing tunneling rate J can be roughly divided into 
two stages, see Fig. [1] The first stage is non-adiabatic 
or impulse in the sense that the state of the system does 
not change and remains close to the uncorrelated initial 
state. After the instantaneous transition rate falls be- 
low the instantaneous relaxation rate of the system, at a 
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J, the evolution crosses over to the adiabatic stage. In 
this stage the phases become increasingly correlated. Af- 
ter the impulsc-adiabatic crossover at J the system is in 
(local) thermal equilibrium. Both its temperature and 
excitation energy scale like the third root of the transi- 
tion rate, apparently in agreement with the experiment in 
Ref. Q. However, in 2 and 3 dimensions even in the adi- 
abatic stage the system has not enough time to develop 
the (quasi-)long-range order expected in thermal equilib- 
rium at low temperature. This lack of long-range order 
manifests itself in the most spectacular way by topologi- 
cal vortex excitations. 

The paper is organized as follows. In Section |TT] we 
introduce the Bose-Hubbard model and define the linear 
quench (ramp) of the tunneling rate from the initial Mott 
state to the superfluid phase. In Section JTiTl we briefly de- 
scribe the truncated Wigner method [15| where bosonic 
operators are replaced by c- numbers, but expectation val- 
ues are obtained as averages over stochastic realizations 
of the initial state. In Section HVl we focus on the Joseph- 
son regime of relatively weak tunneling rate and use the 
truncated Wigner method to investigate thermalization 
of random initial states. We find that there is robust local 
thermalization in all relevant cases and estimate thermal- 
ization time. This is in agreement with previous studies 
of local thermalization after sudden change of tunneling 
parameter [20|. Motivated by this result, in Section fVl 
we consider an adiabatic process in the Josephson regime 
driven by a time-dependent tunneling rate and derive its 
corresponding adiabate equation. Given the estimate for 
thermalization time, in Sections I VII and IVIII we localize 
the crossover between the initial non-adiabatic (impulse) 
stage of the Mott-supcrfluid linear quench to the follow- 
ing adiabatic stage and estimate the initial temperature 
in the adiabatic process. The temperature as well as the 
excitation energy in the adiabatic stage are proportional 
to the third root of the quench rate. This power law is in 
agreement with the experiment Q . In Section I Villi we 
investigate one-particle correlation functions in the adia- 
batic stage. In ID the correlation function is exponential 
as it should be in a thermal state. Its correlation length 
grows like the third root of the quench time. In contrast, 
in 2D and 3D we find that the function has thermal- 
ized at short distance, but it did not have enough time 
to develop the (quasi)-long-range order expected in ther- 
mal equilibrium at low temperature. This may be due to 
the topological defects left behind by the nonequilibrium 
transition at densities that are much higher than what 
would be expected from a thermal equilibrium at a cor- 
responding final temperature. Section llXI completes our 
discussion of the homogeneous case, i.e. without a har- 
monic trap, by extending the quench beyond the Joseph- 
son regime deep into the Rabi regime. Finally, in Section 
IXl we repeat our simulations with a harmonic confine- 
ment and show that, as long as the initial Mott cloud 
extends over many lattice sites, the trap potential does 
not alter the third-root scaling of the excitation energy 
with the transition rate. We conclude in Section IXll 



II. BOSE-HUBBARD MODEL 

The model describes spinless cold bosonic atoms in a 
D-dimensional optical cubic lattice [1, [h], EH, [H[ of L D 
sites numbered by a vector s G Z D . Its Hamiltonian 
reads 

H = -J a L a s2 (—alala s a s + V s ala s j .(1) 

(si,s 2 > s V ' 

Here J is the hopping rate between nearest neighbor sites 
and V s is a trap potential. For the uniform case, when 
V s = 0, we assume periodic boundary conditions. In 
our units the interaction strength is 1 /n, where n is the 
average number of atoms per site. In the thermodynamic 
limit there is a quantum phase transition from the Mott 
insulator to superfluid at J cr ~ n~ 2 . 
We drive the system by a linear quench 

J(t) = t/rq , (2) 

starting in the Mott ground state at J = 0, 

\n,n,n,...,n) , (3) 

with the same number of particles 

n > 1 (4) 

at every site. We end either in the Josephson regime 

J < 1 , (5) 

or extend the linear quench to the Rabi regime J> 1. 

III. TRUNCATED WIGNER METHOD 

When n 3> 1 we can replace annihilation operators 
a s by a complex field (j) s , a s ~ ^fn(j) s , normalized as 
l^ s l 2 = L D and evolving with the discrete Gross- 
Pitaevskii equation 

^ = -JV 2 s + (|0 s | 2 -l)^ s , (6) 
see Refs. Here 

D 

V 2 s = ("AH-** - 2 ^s + 0s-e Q ) (7) 

a = l 

is a D-dimensional Laplacian. Truncated Wigner method 
was used to study dynamics of a quantum phase transi- 
tion in Ref. Q . For alternative approaches not using the 
truncated Wigner method see Refs. HEft 

Quantum expectation values are estimated by averages 
over stochastic realizations of <f> s . Each realization has 
different random initial conditions coming from a Wigner 
distribution of the initial state ((3]): 

MO) = ^ (8) 
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with independent random initial phases 6 S (0). 
We consider kinetic and potential energy, 
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(9) 
(10) 



where V Q S = s +e Q — 0s and the overline means average 
over random initial conditions. 



IV. THERM ALIZATION IN JOSEPHSON 
REGIME 
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In this regime density fluctuations are relatively small, 
> s | 2 « 1, and it is convenient to parametrize 



[l + /s 



(11) 



with real / s and 9 S . After elimination of f s <C 1 in Eq. 
(|6]) we obtain Josephson equations 



d 2 e s 
dt 2 



2 J ^ sir 
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(12) 



where the sum runs over sites s' that are nearest neigh- 
bors of s. The initial conditions are random phases 6* s (0) 



and vanishing velocities (0) = equivalent to vanish- 
ing density fluctuations. 

Since the parameter J could be eliminated from Eq. 
(|T2| by introducing a rescaled time variable u = J x ' 2 t, 
the equations have a characteristic time-scale 



J 



-1/2 



(13) 



In particular, if there is relaxation towards thermal equi- 
librium, then t is the thermalization time. 

Thermalization in case of ID is demonstrated in Fig. 
[2j A thermal state in ID has finite correlation length at 
any finite temperature, hence the thermalization time is 
also finite. In contrast, in 3D at low temperature there is 
long-range order. These infinite-range correlations need 
infinite time to develop. Consequently, in 3D short range 
correlations are quick to thermalize, see Fig. [3l but the 
range of long-range correlations grows roughly like the 
square root of time, see Fig. 0] However, the quick local 
equilibration is sufficient to thermalize local observables 
like the energy density. 

In the Josephson regime / s <C 1 in Eq. (fTTjl. Therefore 
the energy (|9ll0p is approximately quadratic 
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FIG. 2: Thermalization in ID (L — 512). In this figure we 
show correlation functions Cr = i(aja s+ fl) = 4>t<f> s+ R after 
an instantaneous quench from the initial Mott state (|3I8|) at 
J = to a final J <JC 1 in the Josephson regime obtained 
with the truncated Wigner method. As predicted in Eqs. 
(|12I13[1 , the three plots for the widely different final J collapse 
in the rescaled time J l ^ 2 t, proving that r ~ J~ 1//2 is indeed 
the characteristic time-scale in the Josephson regime. For 
large J x ^ 2 t the collapsed plots tend to their equilibrium values 
predicted in Eq. (|43[1 in the Appendix, demonstrating that r 
is indeed the thermalization time. 
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FIG. 3: Thermalization in 3D (128 x 128 x 128 lattice). In 
this figure we show short-range correlators Gi, . . . , Ga after an 
instantaneous quench from the initial Mott state (|3l8p at J = 
to the final J — 0.1 in the Josephson regime obtained with 
the truncated Wigner method. These short range correlators 
Cr do thermalize but, unlike in ID, their relaxation time 
clearly increases with R. There is quick local thermalization, 
but it takes longer time to thermalize the system over larger 
scales. 



by the averages 



At temperature T the fields arc distributed with 
cxp[— (£"kin + E pot )/T] and the system is characterized 
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VI. 



IMPULSE- ADIABATIC CROSSOVER 
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FIG. 4: Thermalization in 3D. In this figure we show the cor- 
relation function Cr(u) = ^(a\a s+ R) = <j>* s 4> s +R for several 
values of the rescaled time u — J x ^ 2 t at a fixed J = 0.01. The 
initial state at u — has phase correlations of finite range. It 



is prepared by a linear quench from the initial Mott state (I3I8P 
at J = to J = 0.01 with tq = 512. After this ramp, this 
low energy initial state thermalizes at the fixed J = 0.01 to- 
wards a low temperature thermal state with long-range order. 
While the short range correlations Cr with small R are quick 
to thermalize, see Fig. O the long-range correlations are slow 
to develop the long range order expected at low temperature. 
As we can see in this figure, the range of Cr roughly doubles 
when the time u increases by a factor of 4 i.e. the correlation 
range grows like the square root of time. 



satisfying the equipartition principle, and total energy 
(E) = (E kin ) + (E pot ) = TL D . (17) 



In order to see when the evolution is adiabatic and 
when it is not, we must compare Q the instantaneous 
transition rate 



dJ 

dt 
J 



1 

Jtq 



(21) 



which is huge at early time and small at late time, with 
the instantaneous relaxation rate in Eq. |)13|) 



J 1 / 2 



t 1 ' 2 
172" ■ 



(22) 



which is negligible at early time and large at late time. 
The two rates are comparable near 
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-2/3 



.1/3 



(23) 



(24) 



In a crude impulse-adiabatic approximation, see Fig. [T] 
after J the evolution is adiabatic, but before J it is im- 
pulse in the sense that the state of the system does not 
change despite changing J and the phases 9g remain as 
random as in the initial Mott state. Thus the random ini- 
tial phases survive until the impulse-adiabatic crossover 
near J, when they become initial conditions for the fol- 
lowing adiabatic evolution. 



VII. EXCITATION ENERGY 



V. ADIABATIC EVOLUTION 

When the process driven by the time-dependent J in 
Eq. ([2]) is adiabatic, then the system follows thermal 
equilibrium with a time-dependent T. On one hand, 
due to the changing temperature T(t), its thermal en- 
ergy (E) = TL D changes at the rate 



dt 



(E) 



dT 
~dt 



-L 



D 



(18) 



On the other hand, the same energy changes due to the 
time-dependent Hamiltonian with the time-dependent J 
at the rate 



dt 



(E) 



dJ_ (E kin ) 
dt J 



(19) 



Equating the two rates, (|18p and (|T^|) . we obtain a simple 
equation 



dT 
~dJ 



T 
2J 



(20) 



Therefore T = AJ 1 / 2 , with an integration constant A 
depending on initial conditions, is the adiabatc equation 
describing the adiabatic process. 



The impulse-adiabatic crossover takes place in the 
Joscphson regime when J< 1 or, equivalcntly, for slow 
enough quenches with 

tq » J- 3/2 . (25) 

Since at J the phases remain random, the kinetic energy 
in Eq. (JUJ) is (E kin ) ~ ^^-JL D . Comparing this with 
(E k in) = \TL D in a thermal state, we obtain the initial 
temperature T ~ 47T 3 D J for the adiabatic process begin- 
ning at J. This initial condition determines the constant 
A in the adiabate equation T = AJ 1 ' 2 as A ~ ^^J 1 / 2 . 
Consequently 



T 



jl/2 jl/2 „ jl/2 



-1/3 



(26) 



is the time-dependent temperature in the adiabatic pro- 
cess after J. This solution remains accurate as long as 
J< 1. 

In the adiabatic thermal state after J the kinetic and 
potential energies scale as 

1. 



(-Ekin) — (Epot) — 
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J 



1/2^-1/3 j-D 



(27) 



in consistency with the numerical data in Figure O and 
Table H 
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FIG. 5: The figure shows dependence of the kinetic energy 
density (left column) and potential energy density (right col- 
umn) for ID (the upper row), 2D (the middle row) and 3D lat- 
tice (the bottom row). The lattice size L = 4096, 256, 128 in 
ID, 2D, 3D respectively. For large tq > 1 we observe power 
law behavior consistent with the predicted (E) ~ Tq 1 ^ 3 . The 
best fits to the tails of the energy plots (the solid lines) give 
the exponents listed in Table U The plots also demonstrate 
the equipartition (Skin) ~ {E po t} at J = 0.1 as predicted 
in the Josephson regime where the energy is approximately 
quadratic. 
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TABLE I: The best fits to a in (E Mn/pot ) ~ t q 
tent with a — 1/3. Lattice sizes as in Fig. [5] 



are consis- 



VIII. 



CORRELATIONS 



In ID a thermal correlation function is exponential, see 
Eq. (|40|) in the Appendix, 



C R = -(4%+^) = exp(-i?/0 (28) 
with a time-dependent correlation length 



<» T 



J 1 / 2 

JV2 



J 



1/2 1/3 



(29) 



compare Eq. (|4ip . In the above derivation of £ we used 
Eqs. (j26l |27|) . which are valid for thermal equilibrium in 
the Josephson regime, i. e. when J < 1. At J, when 
the phases arc still as uncorrelated as in the initial Mott 
state, the length is comparable to the lattice constant 1, 
but after J it grows like J 1 / 2 : the system is ordering by 
increasing the range of correlations. At any given J in the 

1/3 

range J< J<1 the correlation length scales like Tq , 



i.e., the system is correlated more for slower quenches. 
This prediction is consistent with the numerical data in 
Fig. E 
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FIG. 6: Correlation functions Cr in ID at J — 0.1 for dif- 
ferent quench times tq and the lattice length L — 4096. The 
functions are exponential, as expected in a thermal state in 
the adiabatic stage of the evolution. Their correlation length 



scales like £ ^ 
predicted 1/3. 



Tn with the best fit a = 0.329 close to the 
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FIG. 7: Long-range tail of the correlation functions Cr in 2D 
(lattice size L = 256) and 3D (lattice size L = 128) at J = 0.1 
for different quench times tq . They are in the adiabatic stage, 
in the sense that local observables have equilibrated, but they 
have not had enough time to develop infinite-range (quasi- 
)long-range order expected in the low temperature phase. In- 
stead the correlations have finite range limited by a finite rate 
at which they can spread across the system. The range grows 

1/2 

with tq faster than Tq . 

In ID in thermal equilibrium there is finite correlation 
length £ and, consequently, finite relaxation time. The 
system can reach thermal equilibrium in finite time be- 
cause it needs to order only up to the finite distance £. 
In contrast, in 2D and 3D the correlation function Cr in 
low temperature thermal equilibrium either decays with 
a power of the distance R (quasi-long-range order in 2D) 
or tends to a constant (long-range order in 3D). In either 
case the equilibrium correlations have infinite range. For 
a system initialized with random phases it is impossi- 
ble to build up such infinite-range correlations in a finite 
time proportional to tq . Thus the system does not reach 
thermal equilibrium at all length scales: it is correlated 
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as in a thermal state up to a finite range, but it remains 
uncorrelated at longer distances. The short range ther- 
mal correlations explain the — | scaling of the excitation 
energy, because the energy is a local observable not sensi- 
tive to the long range correlations. This is the main result 
of our paper that is experimentally relevant and strongly 
supported by the numerics summarized in Table I. 

The long range correlations in 2D/3D are presented in 
Fig. [7J These correlation functions show how the system 
is equilibrating in time. We do not have analytic pre- 
dictions for these non-equilibrium functions. Apparently 
they are not exponential, so it would be pointless to fit 
them with an exponent to find correlation lengths, but 
one can roughly estimate that their range defined as, say, 

1 /2 

the R where Cr falls below 0.25, grows faster than Tq . 

What is more interesting, the absence of the long-range 
order or, equivalcntly, finite correlation range leaves open 
the possibility of topological vortex excitations, see Fig- 
ures [8] and [9] at J = 0.1 and J = 1 respectively. In the 
Josephson regime at J = 0.1 the healing length ~ J 1 / 2 
is less than the lattice spacing 1. Consequently, there 
are only small and weakly correlated density fluctuations 
around the average density \4> s \ 2 = 1 and there are no 
dips in atomic density associated with the topological 
vortices in Fig. [8j i.e., their cores are thinner than the 
lattice spacing. In contrast, at J = 1 when the heal- 
ing length becomes longer than the lattice spacing topo- 
logical vortices develop empty cores, compare the up- 
per and bottom panels in Fig. [9l The isolated topolog- 
ical vortices in 2D shown in Figs. [8] and [9] are excita- 
tions above the thermal equilibrium in the Berezinski- 
Kosterlitz-Thouless phase [17j. A different perspective 
on vortex formation is provided by Ref. fl8j . where vor- 
ticity in both two and three dimensional XY model is 
discussed. 



IX. INTO THE RABI REGIME 

The linear quench can be extended beyond the Joseph- 
son regime into the Rabi regime. When 

J > 1 (30) 

the hopping term dominates over the nonlinear interac- 
tion in Eq. ([6|), but the nonlinearity is essential to keep 
the system thermalized. In the Rabi regime (E) w (Ekin) 
Eq. (fT9"|) becomes 

d, dJ (E kin ) 

d- t {Ekin) A—r ■ (31) 

Consequently (-Ekin) J ■ The proportionality factor is 
fixed by the initial kinetic energy right after the Joseph- 
son/ Rabi crossover near J ~ I. This energy is roughly 
equal to the final energy when the system is leaving the 
Josephson regime: (E) ~ t q 1/3 L d , see Eq. (27]) with 
J ~ 1. Thus the dominant kinetic energy in the Rabi 




FIG. 8: Quench on a 2D periodic 256 x 256 lattice in the 
Josephson regime at J = 0.1. Top panel: phase 9 S for a 
quench time tq = 1638.4. Middle panel: phase S for a 
quench time tq — 52428.8. Bottom panel: atomic density 
|</> s | 2 associated with the phase in the middle panel. The top 
and middle panel demonstrate that the range of phase corre- 
lations increases with increasing tq. Average size of domains 
of constant phase is consistent with the range of 2D corre- 
lations in Fig. [7J In the faster quench (top panel), where 
the range of correlations is comparable to the lattice spacing, 
there is plenty of random vorticity. In the slower quench (mid- 
dle panel) , where the range of correlations is much longer than 
the lattice spacing, it is possible to identify smooth topolog- 
ical vortices (marked by the arrows). The associated density 
fluctuations in the bottom panel are small, within 10% of 
the average density \(f) s \ 2 = 1, and very weakly correlated be- 
tween different sites. This is not really surprising because in 
the Josephson regime the healing length ~ J 1 ^ 2 is less than 
the lattice spacing and, consequently, densities at different 
sites are only weakly coupled and vortex cores are less than 
the lattice spacing. 
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FIG. 9: A quench with tq — 3276.8 at J = 1 on a part of a 
2D periodic 256 x 256 lattice. In the top panel the phase d a 
and in the bottom panel the density distribution \4>s\ 2 - The 
three arrows mark an isolated vortex and a vortex-antivortex 
pair. Unlike in the Josephson regime at J = 0.1 where density 
fluctuations are small and vortex cores are less than the lattice 
spacing, see Fig. |8j here at J = 1 the topological vortices 
in the upper panel are associated with clear dips in atomic 
density marked in the bottom panel. 



regime is 



(E) a (£ kin ) 



Jr, 



-1/3 L D 



(32) 



J like 



(-Ekiii) 
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(-E kin )| J 



JL 



D 



(33) 



i.e., the excitation energy does not depend on the quench 
time tq when tq -C 1. This is consistent with the nu- 
merical data in Fig. [5] 



X. LINEAR QUENCH IN A 3D HARMONIC 
TRAP 



In a harmonic trap the discrete Gross-Pitaevskii equa- 
tion © becomes 

^ = -JV 2 s + ^ 2 s 2 ^ s + |0 s | 2 s . (34) 

The initial state at J — has random phases 6* s (0), as in 
the uniform case, and a Thomas-Fermi density profile 



.(o)r = — (^tf 



2 s 2 ) 



(35) 



for sites s inside a sphere of radius Rtf and zero oth- 
erwise. In order to make comparisons easier, we set 



2 

M 



here to have 



1 in the center of 



the trap just as in our uniform calculations. Numerical 
results, collected in Figs. [10] and [Til demonstrate that 
in the trap the excitation energy also scales with an ex- 
ponent close to — | just as in the uniform case. At the 
larger J = 1 there are large density fluctuations, see Figs. 
[TT1and[T2T and plenty of random vorticity, see Fig. [121 
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It is linear in the time-dependent J and scales like Tq , 
compare Fig. [5] The scaling ([32]) is valid provided that 
the impulse-adiabatic crossover takes place in the Joseph- 
son regime: J< 1 or, equivalently, rg > 1. 

For faster quenches with tq <C 1 the impulse stage 
extends into the Rabi regime, where the thermalization 
rate t -1 ~ 1 is set by the strength of the nonlinearity 
in Eq. ([6]) . It becomes comparable to the transition rate 
J~ 1 dJ/dt = l/t at i ~ 1 or, equivalently, J ~ Tq 1 ^> 

1 when the impulse stage terminates. At J the phases 
remain as random as in the initial state and the kinetic 
energy is (-Ekin) I j — JL D . In the following adiabatic 
stage, this dominant kinetic energy scales with increasing 



10 1 100 10000 

T Q 

FIG. 10: Excitation energy above the discrete Gross- 
Pitaevskii ground state at a given J, (E) — Egs(J), as a 
function of tq for a 3D lattice in a harmonic trap. The tails 
(E) — Egs ~ Tg a for large tq ^ 1 are best fitted with the 
exponents: a — 0.31 for J — 0.1, a — 0.32 for J = 1.0, and 
a = 0.31 for J = 10.0. Here the lattice size is 64 3 and the 
initial Thomas-Fermi radius 7?tf = 10. 

The key difference with respect to the uniform case is 
that turning on the tunneling rate J makes the trapped 
cloud expand with respect to the original Thomas-Fermi 



- J = 


0.0 


- J = 


1.0 


- J = 


10.0 




FIG. 11: A cross-section of \<f> s \ along the line s y = s z = 
in a single realization of the statistical ensemble for tq = 10. 
The initial Thomas- Fermi profile at J = spreads into a near- 
Gaussian wave-packet at J = 10. The simulation performed 
for a 3D lattice of 64 x 64 x 64 sites. 



profile see Fig. [IT] In an attempt to isolate the ef- 
fect of expansion from the (uniform) Kibble-Zurek mech- 
anism, we rerun our simulations with a constant initial 
phase 6> s (0) = across the system instead of the usual 
random initial phases <j8j) . The resulting excitation en- 
ergies are shown in Fig. 1131 They are not only a factor 
of 10 2 lower than in the corresponding Fig. [TU1 but also 
their decay with tq is steeper, i.e., with an exponent 
closer to —1/2 than to the "random" exponent — |. This 
is not too surprising, as the kinetic energy density of the 
random phases, ~ J, far outweights average density of ki- 
netic energy in the Thomas-Fermi profile with a constant 
phase, ~ JR^p, for any reasonable Rtf ^> 1. 

The exponent —1/2 for the constant initial phase can 
be explained by an impulse-adiabatic argument again. 
When J « 1 the Thomas-Fermi profile (|35|) is a good 
approximation to the ground state of the discrete Gross- 
Pitaevskii equation. When Rtf 3> 1 its lowest Bogoli- 
ubov excitation is 

( \ 

6<p s (t) = ai sin(wii + y>) jo(^T^) ( 36 ) 

with real amplitude a and frequency io\ = Jn 2 R^p. Here 
jo (x) is the spherical Bessel function. It is a breath- 
ing mode describing radial flows of particles. In a linear 
quench of the tunneling rate, J(t) = t/rQ, the evolution 
is impulse as long as the transition rate dJ ^ dt — \jt is 

— 1/2 

much less than the uj\, i.e., up to J ~ RtfTq . At J 
the wavefunction is still the initial Thomas-Fermi profile 
(|35p with a constant phase, but the profile is no longer 
the ground state of the discrete Gross-Pitaevskii equa- 
tion and its excitation energy with respect to the ground 
state is 




E-E, 



GS 



J Rtf 



R 2 r' 1/2 
n TF T Q 



(37) 



FIG. 12: Linear quench in the 3D trap with tq = 204.8 and 
the initial Thomas-Fermi radius Rtf = 25 at the tunneling 
rate J = 1 for 128 x 128 x 128 lattice. Top panel: phase in 
the x — y cross-section across the center of the trap. Middle 
panel: density in the x — y cross-section across the center of 
the trap. Bottom panel: cumulative column density in the 
direction perpendicular to the x — y plane. The arrows in the 
top and middle panel point to a vortex. 



It decays with tq with the exponent —1/2. It is also 
much less than the corresponding excitation energy at 



9 



the J for random initial phases, (E) — Eqs — JRtfi f° r 
any reasonable -Rtf ^ 1- 





— J = 0.1 

— J = 1.0 

— J = 10 











10 100 1000 



FIG. 13: The same excitation energy as in Fig. [10] but with 
constant initial phases 8 S = 0. For large tq the energy decays 
approximately like Tq ' 5 at J = 0.1 and Tq 0A at J = 10 i.e. 
faster than for the random initial phases in Fig. 1101 The size 
of the lattice is 64 x 64 x 64. 



surprise: the ramp is too fast, or the system too slow, for 
the initial Mott state to adjust to the increasing tunneling 
rate. In this impulse stage phases at different lattice sites 
remain as uncorrelated as in the initial Mott state, but 
the frozen Mott state gradually deviates from the instan- 
taneous ground state. At some point, however, reactions 
of the system become fast enough to catch up with the 
ramp and the non-integrable system thcrmalizcs locally. 
In the following adiabatic process, its excitation energy 

1/3 

(or temperature) scales like Tq . This mechanism is 
quite insensitive to the trapping potential because the 
kinetic energy accumulated in the initial random phases 
typically far exceeds the kinetic energy due to localization 
by harmonic confinement. The absence of thermalization 
at large scale manifests itself by topological vortex exci- 
tations. 



XI. CONCLUSION 

Our results justify the following simple picture. A lin- 
ear ramp of the tunneling rate at first takes the system by 
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Appendix: Thermalization in the Josephson regime 
after a sudden quench in ID 

The Josephson equations (fl~2|) follow from a Hamilto- 
nian 



(38) 



where p s = 9 S . A thermal state is given by a fac- 
torizable Boltzmann probability distribution f(p, 9) = 
f P (p)fg{0) oc exp(—Hj/T). Here f p is a Gaussian and, 
in case of a ID chain, 



fe(9) oc exp 



2J^cos(fl s+ i -6 S )/T 



}exp[2JcosA0 a /T] 



(39) 



For a chain much longer than a correlation length, fg 
can be approximately factorized into a product of dis- 
tributions for independent random phase steps A9 S = 



's+l 



The product can be used e.g. to calculate thermal 
correlation functions 



Cr 



s + R e - 



R-l 

s=0 



iA6, 



h(2J/T) 



Io(2J/T) 



(40) 



where I m is the modified Bcsscl function. The correlation 
length is 



£ = 1/log 



WJ/T) 



4J 
T 



(41) 



Ji(2J/T)^ 

with the last approximation for small T -C 4 J. 

We are in a position now to analyze the instanta- 
neous quench in Figure [2] from J = to a finite J <C 
1. The Mott state is the initial state right after the 
quench. It is characterized by p s (0) = and random 
9 S (0) and, consequently, its average energy per site is 
E- 1Tl /L = 2 J. This energy is conserved in the following 
evolution with Josephson equations as the system ther- 
malizes to a temperature T with average energy per site 
E T /L = ~ + 2 J [1 - C\\. Since E in = E T the tempera- 
ture satisfies 



T 

2J =2 



2J 



1 - 



h(2J/T) 
WJ/T) 



(42) 



Its solution is x - 
relation length £ 
Cr = x~ R : 



4J/T = 2.1312 leading to the cor- 
= l/log(x) = 1.321 and correlators 



Ci = 0.469, C 2 = 0.220, C 3 = 0.103, C 4 = 0.049 , (43) 

compare with Figure [^1 Notice that the asymptotic ther- 
mal correlators after the sudden quench do not depend 
on J«l. 



